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We have developed a thorough and accurate method of determining anharmonic free energies, 
the temperature dependent effective potential technique (TDEP). It is based on ab initio molecular 
dynamics followed by a mapping onto a model Hamiltonian that describes the lattice dynamics. The 
formalism and the numerical aspects of the technique are described in details. A number of practical 
examples are given, and results are presented, which confirm the usefulness of TDEP within ab initio 
and classical molecular dynamics frameworks. In particular, we examine from first-principles the 
behavior of force constants upon the dynamical stabilization of body centered phase of Zr, and show 
that they become more localized. We also calculate phase diagram for *He modeled with the Aziz et 
al. potential and obtain results which are in favorable agreement both with respect to experiment 
and established techniques. 



I. INTRODUCTION 

One of the common goals for first principles calcu- 
lations is the comparison of energies, such as configu- 
rational energies, surface energies, mixing enthalpies or 
lattice stabilities. Usually the comparison is limited to 
total energies. This is appropriate when the effects of 
temperature can be neglected. However, for many prob- 
lems within physics, materials and Earth sciences they 
are not negligible and Gibbs free energy represents the 
proper thermodynamic potential when the temperature 
and pressure are external parameters. The quasihar- 
monic approximation can bridge the temperature gap, 
but there are cases where it falls short. Strongly an- 
harmonic systems are not described well^, especially dy- 
namically unstable systems. Traditionally the problem 
of dynamical instability was addressed either by includ- 
ing more terms in the Taylor expansion of the potential 
energy or via a self-consistent approach .I^Hll 

HootoiP realised that even though the second deriva- 
tives at the equilibrium positions are negative, the atoms 
in a solid rarely occupy these positions. They move in 
the effective potential of their non-stationary neighbours. 
By sampling the potential energy surface not at the equi- 
librium positions but at the most probable positions for 
a given temperature one can obtain a harmonic approx- 
imation that describes the system at elevated tempera- 
tures. The self-consistent formalism employs an iterative 
procedure- by creating a harmonic potential, which is 
used to describe the thermal excitations that again give 
a new harmonic potential. This is then repeated until 
self-consistency. 

The double-time Green's functions, developed by 
Choquard,^ use a cumulant expansion in the higher order 
terms. Although formally exact, this formalism requires 
knowledge of the higher order force constants. Obtain- 
ing these accurately for something other than the struc- 
turally simple systems is computationally very demand- 
ing from first principles.'^ 

A recent implementation of the self-consistent formal- 



ism by Souvatzis et al^^^ uses ah initio supercell cal- 
culations. A problem with this approach is that the 
excitations could only be done in the harmonic sense 
which means probing phase space with a limited basis set. 
Where the harmonic approximation works well this is not 
a problem. When the harmonic approximation fails it is 
due to a strong anharmonic contribution. Strongly an- 
harmonic systems are by definition badly described with 
the harmonic approximation. 

Several techniques use Born-Oppenheimer molecular 
dynamics to obtain anharmo nic corrections to quasi- 
harmonic free energiesl ^^ l ^^ l ^^ ^They focus on anharmonic 
corrections to materials that are well described in the 
quasiharmonic approximation, and the applicability to 
strongly anharmonic systems is questionable when the 
phonon renormalization due to increased temperature 
can not be described by a linear equation. 

We have developed a new method^ that is similar to 
Hootons^ original idea, but with a foundation in (ab ini- 
tio) molecular dynamics. In this paper we present a sub- 
stantial refinement and generalization of the temperature 
dependent effective potential method (TDEP), showing 
how it deals with: a model one-dimensional anharmonic 
oscillator, a strongly anharmonic system, bcc Zr, treated 
from first-principles, and ^He modelled with the Aziz et 
al. potential. 

II. TDEP FORMALISM 

The starting point of our method is to introduce a 
model Hamiltonian for a Bravais lattice in the harmonic 
form: 

i * ij 

which describes the lattice dynamics. Here and are 
the momentum and displacement of atom j, bold symbols 
indicate vectors and doubly overlined symbols matrices 
respectively. The reference point for the displacements 
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is the OK relaxed lattice (initially, in Sec. IV we will 
revisit this) . The interatomic force constants § and the 
ground state energy Uq are yet to be determined. Given 
Na atoms in this model system the forces acting on the 
atoms are given by 
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As detailed in our previous papeJ^^ we seek to determine 
the force constant matrices through minimization of the 
difference in forces from the model system a real system, 
simulated, for instance, by means of ab initio molecular 
dynamics (AIMD). An AIMD simulation will provide a 
set of displacements {Uj^^}, forces {F^^} and potential 
energies {Ef^^}. We seek to minimize the difference in 
forces from AIMD and our harmonic form (F^) at time 
step t, summed over all time steps Nt: 
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This is realized with a a Moore-Penrose pseudoinverse 
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to obtain the linear least squares solution for $ that min- 
imize AF. This is then mapped to the form 
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where a/3 are indices to atoms in a unit cell with Nuc < 
Na atoms and fj-v Cartesian indices. The pair vectors 
in the supercell TLij are mapped to stars of lattice vec- 
tors R; connecting atoms of type a and /3. From this 
form the phonon dispersion relations, free energy and all 
other quantities can be extracted. This direct implemen- 
tation works well,li3 but the numerical efficiency can be 
improved, as is demonstrated below. 



III. SYMMETRY CONSTRAINED FORCE 
CONSTANT EXTRACTION 

The form of the force constant matrices depends only 
on the supercell used in the AIMD and the crystal lat- 
tice. We begin by populating the force constant matrices 
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and so on, including vectors R; up to a cutoff deter- 
mined by the size of the supercell. Some of these 9k are 
equivalent. To find out which, we look at the symmetry 
relations the force constant matrices have to obey. We 
havcP 
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that stem, in order, from the facts that there is no net 
translation of the crystal, all Bravais lattices have inver- 
sion symmetry, and that the order of the second deriva- 
tives does not matter. Each relation will give us a few 
equations for the unknowns 9k, reducing their number. 
In addition to these fundamental properties of the force 
constant matrices we can benefit from the symmetry of 
the lattice. If symmetry relation S, belonging to the 
point group of the lattice, relates two vectors Rj = 5'Rfc 
we have the following relation: 



* ^(R0 = 5* (Rfc)^' 



(10) 



By applying equations ([7l)-(|10p the number of unknown 
variables is massively reduced. For example, a bcc lat- 
tice modelled as a 4 x 4 x 4 supercell (128 atoms) would 
have 147456 unknown variables in if one docs not con- 
sider symmetry arguments. Application of Eqs. (|7|-( |10[ ) 
reduces the problem to 11 unknown variables. Having 
found the reduced problem with Nq unknown variables, 
it can be substituted back into ([2|. The expression for 
the forces at timestep t will schematically look like this: 
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The actual distribution of the 9k will depend on the prob- 
lem at hand. Carrying out the matrix product gives us 
a new set of equations for the forces 



(12) 
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where second sum describes the coefRcients for each 
9k contained in the expression for force component 
These coefRcients are hnear combinations of the displace- 
ment components us- The exphcit form is determined by 
the lattice. In matrix form this is written as 



F = c(u)©, c(u)fc^ = ^c: 



(13) 



Eq. 13 is equivalent to Eq. [2] it is just rewritten in terms 
of the symmetry inequivalent interactions. This imple- 
mentation symbolically reduces the number of unknowns, 
generates the function that gives the matrix C from a set 
of displacements U and the mapping from the set of 9^. 
back to the force constant matrix <^'^^{Tli)^^ 

To find Q we minimize the difference in forces between 
AIMD simulations and the model hamiltonian 
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This is again realised with the least squares solution for 
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Having determined we can substitute the components 
into the force constant matrix and proceed to calculate 
thermodynamic properties of the original (real) system. 

The suggested scheme is a superior way of using sym- 
metry to improve the numerical accuracy. Most schemes 
involving symmetry revolve around determining the in- 
teraction in one direction and then using the symmetry 
relations to translate and rotate that interaction. Any 
numerical errors-which are always present-will propa- 
gate to all interactions, whereas in the present approach 
the errors will be averaged and should cancel each other 
to a large degree. 

The advantage of this procedure will be illustrated be- 
low. 



IV. INTERNAL DEGREES OF FREEDOM 

If the system has internal degrees of freedom for the 
structural relaxations, such as a crystal with point de- 
fects, an interface or a random alloy the atoms ideal po- 
sitions and equilibrium positions do not coincide. While 
one could find the relaxed positions from OK calculations. 
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FIG. 1. (color online) Convergence of the free energy of bcc 
Zr at 1300K with respect to the number of timesteps. In panel 
(a) symmetry is treated numerically and in (b) it is treated 
in the novel analytical way. The same input data is used in 
both cases, and it converges to the same value, but we have 
improved the performance by several orders of magnitude. 



the equilibrium positions are by no means constant with 
respect to temperature. TDEP handles this in an el- 
egant way. Note that we find the second order terms 
in ([T]) with a least squares fit of the slope of force versus 
displacement. Originally, the displacements could be cal- 
culated with respect to the wrong equilibrium positions 
that do not correspond to the temperature of the simula- 
tion. Still our experience shows that slope will be the well 
approximated. That allows for the following procedure: 

- Guess equilibrium positions, usually the ideal lat- 
tice positions. 

- Use these to calculate the displacements u from the 
AIMD simulations. 



Determine and from that the residual force 

Nt 
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where Fj^^ are the AIMD forces and F^ are given 
by equation |13[ N^ is the number of timesteps, and 
subscript t denote the corresponding forces at time 
step t. 

- These forces are then used to move the atoms in a 
steepest descent scheme towards equilibrium posi- 
tions. The whole process is repeated until conver- 
gence. 

Our test shows that this procedure is remarkably sta- 
ble. The second order force constants $ are weakly af- 
fected by the choice of equilibrium positions. The vibra- 
tional entropy and phonon dispersion relations are largely 
unaffected as well. Eliminating the first order term, how- 
ever, is formally important and crucial when extract- 
ing higher order terms. Note that self-consistent itera- 
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tions are numerically efficient, because the most time- 
consuming step for applications of TDEP is MD simu- 
lations, while their mapping on model Hamiltonian ([T]) 
represents post-processing of the MD results with mini- 
mal computaional cost. 



V. DETERMINING THE FREE ENERGY 

We will begin by reiterating the traditional way how 
free energy is determined in the quasiharmonic approxi- 
mation. Divided into parts it will be 
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where a division of the right hand side into parts makes 
a clear distinction between the electronic contribution 
Fci and the vibrational contribution i^vib- The electronic 
contribution is divided into the total energy of the lattice, 
Utot , and the electronic entropy Sei ■ The vibrational con- 
tribution is divided into average kinetic Ek and potential 
energy t/vib of the ions, vibrational entropy S'vib and zero 
point energy U^p. The lattice contribution is obtained 
from DFT calculations with the Mermin functional and 
the vibrational part from the harmonic approximation 
via 
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Where giuj) is the phonon density of states. In this ap- 
proach, all the vibrational contributions are calculated 
within the harmonic approximation. 

Turning to AIMD, the free energy (in the canonical 
ensemble) is divided as 



{Umd) + {Ek} — TSmd, 



(19) 



were the potential energy C/md is temperature dependent. 
Since the ions move as classical particles the zero point 
energy is missing. There is unfortunately no informa- 
tion about the entropy, but through the force constant 
matrices obtained using TDEP, the vibrational entropy 
and zero point energy can be estimated. For TDEP to 
have an accurate free energy the potential energy should, 
on average, be equal to that of AIMD. This would ensure 
that the full anharmonic (Umb) is included. The problem 
is that Umd is rapidly oscillating over time, requiring a 
long simulation time to converge. If we look at the TDEP 
potential energy 
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and recognise that it should model the thermal fluctua- 
tions of the original system we can overcome the numer- 
ical issues. Setting the average potential energies equal. 
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FIG. 2. Illustration of the different terms included in free en- 
ergy calculations using different approaches. The filled boxes 
denote the terms that are included, and the striped areas what 
is omitted. The main message is to point out that the internal 
energy has a non-trivial temperature dependence, something 
that is omitted in the quasiharmonic approximation. HOT 
indicate the higher order terms that are missing in the TDEP 
free energy, Eq. [22] and [18] 
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By removing the thermal excitations of the harmonic 
form the fluctuations can be decreased by roughly an or- 
der of magnitude and the accuracy of Uq is thus increased 
be the same amount. Including higher order terms in the 
energy expansion would further serve to minimize these 
fluctuations. 

The potential energy that was removed will be added 
again when the Helmholtz free energy is calculated: 
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where i^vib is the phonon contribution given Eq. [18] with 
the phonon density of states determined in the TDEP 
formalism. It includes the kinetic and potential energy 
of the ions. In Fig. [2] we further illustrate the difference 
in methods of obtaining the free energy. 

The formally exact method of thermodynamic 
integratiorf^ can be used to determine the free energy. 
This method will determine the anharmonic correction 
to the free energy. If the TDEP model Hamiltonian is 
used as the reference point the full free energy will be 

F^Ua+ F,ib + / {Umd - (7tdep)a d\ . (23) 
Jo 

" V ' 

AFah 

The integral is over the Kirkwood coupling parameter A, 
and the potential energy difference is between the model 
Hamiltonian and the molecular dynamics potential en- 
ergy. 

The model TDEP Hamiltonian is constructed to de- 
scribe the system as accurately as possible while retain- 
ing the harmonic form. It is then easy to argue that 
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the anharmonic correction term A^ah should be small. 
While it is difficult to make general statements regarding 
this, our experience so far is that this correction is very 
small in every system we have tested. 

In addition, the thermodynamic integration technique 
can be numerically inefficient when high accuracy is 
needed. While one can accurately control the numerical 
accurac jffSlthe finite size effects are more difficult to con- 
trol, especially in ab intio simulations. In Fig. [3]we show 
that the error due to the limited size is on the same order 
as the correction to the TDEP free energy in reasonable 
simulation size^^^. It makes little use to add a correction 
to the TDEP free energy where the uncertainty is of the 
same order of magnitude as the correction itself. 

The TDEP free energy, on the other hand, behaves 
well with respect to limited simulation cell. In figure |4] 
we see that at the reasonable system size of 100 atoms the 
free energy is converged within ImeV/atom.'^^l It is also 
easily converged in terms of simulation length: in Fig. [5] 
we illustrate the advantage of analytically treating sym- 
metry. In our previous work,'^^ we studied Zr in the bcc 
phase. There, we found convergence within ImeV/atom 
for the free energies after 25000 time steps. With the 
symmetry constraints, we converged to the same value 
using 50 time steps, an improvement by several orders of 
magnitude. 
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FIG. 4. (color online) Convergence of the free energy from 
TDEP with respect to simulation size. The system in question 
is fee Cu with a classical embedded atom potential^. At sizes 
about 100 atoms it is converged within ImeV/atom. 
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FIG. 5. (color online) Convergence of the free energy of bcc 
Zr at 1300K with respect to the number of timesteps. In panel 
(a) symmetry is treated numerically and in (b) it is treated 
in the novel analytical way. The same input data is used in 
both cases, and it converges to the same value, but we have 
improved the performance by several orders of magnitude. 



FIG. 3. (color online) Convergence of the free energ y co rrec- 
tion from thermodynamic integration ( A_Fah in Eq. 23 1 with 
respect to system size. At system sizes smaller than ~700 
atoms the uncertainty is of about the same order as the cor- 
rection. This particular case is for fee Cu modelled with an 
embedded atom potentiaP^. 



VI. APPLICATION OF TDEP TO A 
ONE-DIMENSIONAL ANHARMONIC 
OSCILLATOR 

To illustrate TDEP we first apply it to a one- 
dimensional anharmonic oscillator. Consider the follow- 
ing one-dimensional potential: 
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Here x is the position and fc, a and /? are known param- 
eters. The equilibrium position can depend on temper- 
ature and is assumed to be at T = OK. The aim is 
to find the second degree polynomial fit to Eq. [24] that 
best describes the system. If this polynomial only con- 
sist of a quadratic and a constant term it will describe a 
harmonic oscillator with well defined free energy. If one 
applies the harmonic approximation to this potential it 
will not work well. The second derivative 



dx^ 



= k~ 2a/3e~^(^-^")' + (4a/3(a; - xo)fe-^^''-''°^' 



(25) 

will determine the force constant $. The temperature 
dependence of xq is omitted and we will end up with 
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This, as seen in Fig. [6j will not be a particularly good 
model for the true potential. These issues arise from the 
fact that the potential energy surface is only probed at 
a; = 0, the T = equilibrium positions. To work around 
this problem, let us apply TDEP and put a particle in 
the potential given by equation ( 24 ) and perform a MD 
simulation. When controlled by an appropriate thermo- 
stat, the particle will yield a set of Nt forces, positions 
and energies, {Ft,Xt,Et}, one for each timestep. This 
data can now be used to fit a potential of the form 
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FIG. 7. Comparison of the effects of including higher or- 
der terms between the harmonic approximation and TDEP. 
When extending the Taylor expansion of an anharmonic po- 
tential (dashed blue line) in the Born-von Karman ansatz to 
higher order terms we end up with the series of lines depicted 
in panel (b). Panel (a) shows the same extension for TDEP. 
Even limiting oneself to the second order term {n — 2), the fit 
will implicitly contain anharmonism to an arbitrary degree in 
the range that is thermally accessible. Extending TDEP to a 
higher order converges towards the true potential faster than 
when higher order terms are added to the harmonic approxi- 
mation. 



FIG. 6. Comparison of performance of TDEP and the har- 
monic Taylor expansion for the potential described by equa- 
tion |24[ Three examples are shown when the conventional 
harmonic approximation fails to describe the potential while 
TDEP succeeds. Panels (a),(c) and (e) show the potentials 
and (b), (d) and (f) show the forces, a and xo are the param- 
eters of Eq. 
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In panel (a) the harmonic approximation cor- 
respongd to a dynamically unstable system, whereas TDEP 
provides a dynamically stable solution. In panel (b) the har- 
monic approximation provides an inaccurate potential. Panel 
(c) shows how TDEP finds the high temperature equilibrium 
position Xo- 



U{x] 
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a harmonic potential centered at x'q. Let us begin by 
determining the second order term. As discussed in Sec. 
IV we guess a value for x'q , use the forces from molecular 
dynamics {-Ft} and minimize 
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This is easiest realized as a least squares fit of a straight 
line in forces, as demonstrated in the right panels of 
Fig. [6j Equation ( 28 1 determines the second order term. 
The residual force at Xq, AF, can be used to find the 



equilibrium position. It is done in the following manner: 



a guess for x'q gives us a and A-F. This residual force 
is used to move x'q to a new position, and the process is 
repeated until self-consistency is reached. When we have 
found the equilibrium position we can safely assume that 
any first order term in our polynomial can be set to 0. 
As described in Sec. |v]The constant energy term, <^^'^\ 
can be determined from the energies {Et \ obtained from 
molecular dynamics simulations: 



$(0) = 



(E,-\^^-\x,-x',fY 



(29) 



This is the best possible potential of the harmonic form 
at a given temperature that approximates the original 
potential in Eq. [24] In Fig. [6] the true potential and 
the fit are illustrated for different a,j3 and xq. The an- 
harmonism of the potential is implicitly described by the 
polynomial fit. In Fig.[7]the expansion in Eq. [27]has been 
extended to higher orders for an anharmonic potential. 
TDEP, probing the effective potential at finite temper- 
ature, converges to the true potential rapidly whereas 
including more terms in the Taylor expansion in Eq. ([T]) 
does by no means guarantee numerical stability at finite 
temperature. 
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VII. PRACTICAL APPLICATIONS OF TDEP 

Let us summarize this and present the scheme used 
to calculate accurate Gibbs free energy surfaces in the 
TDEP formalism from first principles. 

- First calculate DFT total energy as a function of 
volume. This provides an equation of state and 
allows us to choose the volume interval, which 
covers pressure of interest. 

- If feasible, obtain approximate harmonic 
potentials for the systems at hand in this volume 
interval. These potentials are used to speed up 
the calculations as described in Steneteg et at 
(submitted to PrB, BZ11707). 

- On the grid of volumes and temperatures, perform 
AIMD simulations in the canonical ensemble. 

- From these simulations, extract internal energy Uq 
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(21 1 and interatomic force constants using Eq. 14 
ensuring convergence of the free energy with 
respect to simulation length. 

- To increase further the accuracy of the 
calculation, it is recommended to select a subset of 
uncorrelated samples from the AIMD simulations 
upsample these to high accuracy, as described in 
Ref. [m A new free energy is calculated. 

- The equation of state is interpolated over the grid 
of temperatures and volumes providing the Gibbs 
free energy surface. This is then repeated for each 
structure, compound or composition of interest. 

TDEP is a thorough and time-consuming method, but 
the results are excellent. The phonon dispersion rela- 
tions of a material that is dynamically unstable at zero 
temperature is a good example. When re-evaluating the 
results for Zr obtained in Ref. 13, we observe a striking 
difference in harmonic and TDEP force constants. This is 
illustrated in Fig. [Sj The effective TDEP force constants 
decrease faster with distance compared to the harmonic 
ones, a behavior that is expected. It is a vivid illustra- 
tion of the temperature dependence of the potential en- 
ergy landscape, and at the same time a confirmation that 
the TDEP technique describes this renormalisation well. 
From the free energy surface we can extract the finite 
temperature equation of state for bcc Zr, as illustrated 
in Fig. [9] 

To test the performance of TDEP close to melting, 
we turn to solid He modelled with the Aziz et al. 
potentia P^ I I^° l The melting curve and bcc-fcc transi- 
tion just before melting has been extensively studied, 
see Ref. and references therein. As demonstrated in 
Fig. [To] strongly anharmonic He pose no problem for the 
presented method. The stabilization of the bcc phase 
before melting is consistent with results from phase- 
coexistence simulations. This are at the moment consid- 
ered the most accurate methods for determining phase 



o 



0.5 



o 











I 


I p I 





2 3 

Coordination shell 



FIG. 8. (color online) Comparison of the TDEP and quasi- 
harmonic force constant matrices. We have plotted the Frobe- 
nius norm of the force constants versus coordination shell. In 
the inset we show the corresponding phonon dispersion re- 
lations. The empty circles are OK harmonic values and the 
filled circles are TDEP values extracted at 1300K. At high 
temperature the interactions at close distances are stronger, 
and fall ofi' faster with increasing distance. 



stabilities at high temperature. They do, however, re- 
quire simulation cells much larger than what is accessible 
to AIMD, and can only be used with classical potentials. 
We show here that with the presented method we can 
with simulation sizes of 125 atoms accurately reproduce 
the same transition temperatures. This verifies the ac- 
curacy of the method and opens up the applicability to 
high pressure high temperature studies of phase stabili- 
ties close to melting. 



VIII. CONCLUSIONS 

We have presented detailed description of the Tem- 
perature Dependent Effective Potential method for the 
treatment of lattice dynamics of strongly anharmonic 
solids, including an extension and refinement to this ac- 
curate technique. Moreover, we have detailed how the 
temperature dependence of all components of the free 
energy should be taken into account, and presented sev- 
eral successful examples, including a model anharmonic 
potential, first-principles calculations of equation of state 
for bcc Zr, and classical molecular-dynamics simulations 
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FIG. 9. (color online) Equation of state for bcc Zr calcu- 
lated with TDEP method at temperatures 500K (dashed line), 
lOOOK dot-dashed line, and 1500K (dotted line). Equation of 
state obtained by means of conventional DFT calculations at 
T=OK is also shown with full line, and the zero temperature 
equilibrium volume Vo=22.83A^ is chosen as a reference point 
at all the temperatures. Note, that bcc Zr is dynamically un- 
stable at T=OK, see the bottom inset in Fig. 7. 



FIG. 10. (color online) The calculated phase diagram for *He 
modelled with the Aziz et al. potential. The red line indi- 
cates the experimental melting curve. The observation of the 
stabilization of the bcc phase before the melting demonstrates 
that TDEP treats this system accurately and in agreement 
with other approaches even in such a strongly anharmonic 
system as ^He. 



of bcc-to-fcc transition in "^He. 
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